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FINITE  DIFFERENCE  MODELLING  OF  SEISMIC  WAVE  PROPAGATION 

AT  THE  SEAFLOOR 


Final  Technical  Report  for  NORDA  Contract  #N000I4-88-M-66]0 


R.A.  Stephen 
RASCON  Associates 
P.O.  Box  567 
W.  Falmouth,  MA  02574 


1.  INTRODUCTION 

The  objectives  of  this  study  were  i)  to  implement  the  WHOI  finite  difference  code 
on  the  CONVEX  computer  at  NORDA  and  ii)  to  run  a  suite  of  models  on  the  effects  of 
lateral  heterogeneity  on  the  primary  response  from  the  seafloor. 

Finite  difference  solutions  to  the  elastic  wave  equation  accurately  predict  the 
response  of  impulsive  and  continuous  wave  sources  in  media  with  arbitrary  vertical  and 
horizontal  variation,  with  fluid-solid  interfaces  ajid  with  shear  propagation  in  the  solid.  All 
possible  wave  types  are  included  (reflections,  refractions,  diffractions,  Stoneley  and 
pseudo-Rayleigh  waves,  evanescent  waves  and  head  waves).  The  primary  disadvantage  of 
the  method  is  that  it  is  very  computation  intensive  and  it  is  generally  limited  to  problems 
with  dimensions  of  only  a  few  tens  of  wavelengths.  This  can  be  partially  alleviated  by 
using  powerful  computers  such  as  the  CONVEX  at  NORDA. 

This  proposal  addressed  the  following  categories  in  NORDA-BAA-88-2; 

1)  Acoustic  AS W  oceanography 

b)  ultra  low  and  very  low  frequency  propagation 

i)  acoustic  transients 

j)  low  frequency  arctic  acoustics 

o)  acoustic  field  interactions  (scattering) 

2)  Non-acoustic  oceanographic  measurements 

i)  bathymetry 

j)  geophysics  i 

5)  Computer  Modelling 

e)  geoacoustic  extensions  to  environmental  acoustic  prediction  systems. 

The  latter  item.  5)e),  is  the  particular  focus  of  the  proposal. 


2.  THE  FINITE  DIH'ERENCE  METHOD 


The  finite  difference  technique  (Stephen,  1988)  is  a  powerful  method  for  studying 
tlie  complete  sea  bottom  interaction  of  the  acoustic  field  in  the  ocean.  Because  llie 
computational  effort  is  quite  large  the  metliod  is  generally  restricted  to  low  frequencies 
(5-25  Hz)  but  in  many  cases  it  is  at  these  frequencies  where  bottom  interaction  becomes 
imporurnt.  Because  the  metliod  is  based  in  the  time  domain  it  is  particularly  well  suited  to 
pulse  or  transient  problems.  The  technique  also  has  considerable  promise  for  studying 
scattering  from  ice  in  arctic  acoustics.  Tie  code  lends  itself  well  to  extending  acoustic 
forward  modelling  schemes  to  include  geoacoustic  bottom  properties. 

Stephen  (1988)  has  reviewed  the  finite  difference  method  as  applied  to  bottom 
interaction  problems,  A  number  of  applications  including  focusing  of  detenninistic 
stnjcture  (Stephen,  1984;  Dougherty  and  Stephen,  1987;  Stephen,  1988)  and  scattering 
from  random  seafloor  (Dougheny  and  Stephen,  1988)  have  been  presented.  Also  the  code 
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(Stephen,  1983)  and  by  participating  in  the  benchmark  sessions  at  the  Acoustical  Society  of 
America  meetings  (Stephen,  submitted).  We  have  sufficient  exp>ericnce  with  the  code  that 
we  are  confident  that  it  wall  be  useful  in  studyiiig  bottom  interaction  problems  in  range 
dependent  environments  including  shear  wave  effects  in  tlie  bottom.  NORDA  should  be 
interested  in  the  results. 


3.  IMPLEMENTATION  OF  THE  FINITE  DITTERENCE  CODE  AT  NORDA 

Tc  objective  here  was  to  install  the  WHOl  finite  difference  code  on  the  CONVEX 
and  to  run  test  models  at  lOHz  out  to  10  km  for  a  full  ocean  depth  of  5.5  km.  Deep  models 
like  this  had  not  been  run  before  so  w'e  made  a  number  of  modifications  before  we  could 
nm  this  model  successfully.  Te  models  are  described  in  section  4.  Here  w'e  review  the 
modifications  tliat  were  made  and  describe  where  the  files  are  on  the  CONVEX. 


A  finite  difference  model  is  run  in  three  steps:  i)  a  preprocessor  which  sets  up  the 
necessary  arrays  for  a  given  calculation;  ii)  a  program  to  compute  the  elastic  parameters  and 
density  in  a  transition  zone  near  the  seafloor,  and  iii)  the  actual  finite  difference  calculations. 
Software  for  each  stage  including  makefiles  for  compilation  is  located  in/mntystephen/prep, 
/mnt/stephen/bny,  /mnt/stephen/diff  respectively  on  the  CONVEX.  A  macro  (or  command 
file)  for  each  model  run  is  located  in  that  model's  subdirectory.  For  example,  the  macro  for 
nuxlel  BBNYl  is  located  in/obs3b/stcphcn/bbnyl/bbnyl.bch.  Once  the  parameter  file  is 
defined  a  nxxiel  can  be  run  by  submitting  (he  .bch  file  to  the  batch  queue.  More  details  on 
the  structure  and  implementation  of  the  code  can  be  found  in  Hunt  et  al  (1983). 

In  /mnt/stephen/diff  are  all  the  versions  of  the  finite  difference  code  used  in  this 
study.  SFINDIF.FOR  is  the  main  driving  code  for  all  cases.  The  changes  occur  in  the 
subroutines  SDUMTS*  and  SBIBTS*,  which  actually  c-irry  out  the  template  calculations, 
and  in  the  subroutines  SFINSUB*  and  SBIBSUB*  which  carry  out  the  absorbing 
boundary'  calculations. 

We  used  two  templates.  Initially  we  used  tlic  Bhasavanija  template  (Stephen  et  al, 
1585;  Stephen,  1988)  however  this  resulted  in  instabilities  at  the  fluicl'solid  interface  at 
large  times.  The  second  template  is  ba.sed  on  a  fomiulation  presented  by  Virieux  (1986). 
lliis  produces  stable  and  accurate  results  for  a  wide  range  of  models. 

We  used  two  absorbing  boundary  schemes.  The  first  w  as  based  on  e  parabolic 
equation  approximation  right  at  (he  boundary  and  the  second  was  based  on  the  telegraph 
equation  applied  in  a  region  near  the  boundary.  The  telegraph  equation  has  temis  which 
introduce  attenuation  into  the  grid. 

The  original  WHOI  finite  difference  code  used  SDUMTS5  and  SF1NSUB4.  Tliis 
used  a  Bhasavanija  tciiiplate  with  parabolic  equation  boundaries  on  the  right  and  bottom 
sides.  (The  top  side  is  a  free  surface  and  the  left  side  is  an  axis  of  symmetry.)  Although 
adequate  for  small,  short  runs  this  was  unstable  on  the  bottom  edge  for  the  large  models. 


The  next  try  (SBIBTS5  and  SBIBSUB4)  used  the  telegraph  equation  on  the  bottom 
edge  only.  This  was  unstable  on  the  right  boundary. 

The  third  case  (SBIBTSG  and  SBIBSUB4)  used  the  telegraph  equation  on  the 
bottom  edge  and  the  right  hand  edge.  The  grid  boundaries  were  stable,  but  an  instability 
occurred  at  the  seafloor  at  large  times.  After  all  of  tlie  principle  phases  had  passed  ilirough 
the  seafloor,  a  phase  app<iared  at  the  interface  which  looked  like  a  Sioneley  wave. 

However  its  amplitude  grew  unrealistically  with  time.  We  call  this  the  'Bibee'  instability. 
We  do  not  recommend  using  SDUMTS5,  SB1BTS5,  SB1BTS6,  SFINSUB4  and 
SBIBSUB4  for  large  models.  They  are  on  the  CONVEX,  however,  for  completeness. 

In  order  to  avoid  the  'Bibce'  instability  we  went  to  a  Virieux  fonnulation  with  the 
telegraph  equation  on  the  bottom  edge  and  parabolic  equation  on  the  right  edge  (SBIBTS7 
and  SBIBSUB7).  This  fixed  the  'Bibce'  instability  but  there  were  false  reflections  from 
the  right  side. 

Finally  we  used  the  telegraph  equarion  on  the  right  and  bottom  sides  (SBIBTS9  and 
SBIBSUB7).  This  works  fine  and  the  test  imodels  in  the  next  section  use  this  fonnulaiion. 

We  also  modified  the  WHOl  code  to  output  run  time  infonnation  and  ma.\imum 
amplitude  values  to  the  log  file  (*.LG4)  during  execution.  Tliis  facilitates  testing  and 
provides  a  record  of  the  run  time  for  each  job. 

While  testing  we  used  the  snapshot  display  on  the  SUN  developed  at  NORDA. 
However,  time  series  are  more  physically  significant.  To  aid  in  the  reading  and  plotting  of 
the  time  scries  files  (*.TST)  we  wrote  a  program,  askii.f,  which  reads  the  binary  file 
U'.TS'F)  and  writes  an  askii  file  (*.ASK).  It  is  possible  in  askii.f  to  siibsample  the  receiver 
locations.  Askii.f  can  be  used  as  the  basis  for  any  code  that  needs  to  read  a  *.TST  file, 
such  as  plotting  code  or  a  seismic  processing  code. 


4.  SUITE  OF  TEST  MODELS 


In  order  to  conrirr!i  that  die  finite  difference  code  on  the  CONVEX  actually  solved 
problems  of  interest  to  NORDA,  we  ran  a  suite  of  six  test  models: 

i)  a  range  indcix^ndent  layered  model  representative  of  the  seafloor  (BBNYl); 

ii)  a  model  like  i)  but  with  bottom  roughness  (BBNY2); 

iii)  a  model  like  i)  but  with  basement  roughness  (BBNY3); 

iv)  a  model  like  i)  but  widi  both  bottom  and  basement  roughness  (BBNY4); 

v)  a  model  like  i)  but  with  a  discontinuous  high  velocity  stringer  in  the  sediment 
(BBNY5); 


vi)  a  model  like  i)  but  with  a  different  shear  velocity  profile  in  the  sediments  and 
basement  (BBNY6). 

Each  model  on  the  CONVEX  has  its  own  subdirector>'  (/obs3b/stcphen/bbny*). 

The  model  is  lun  by  submitting  bbny*.bch  to  the  batch  queue.  This  file  must  be  nnodified 
to  include  the  correct  *.PAR  file  (change  all  occurences  of  BBNY*  to  the  same  correct 
name).  Tlie  makefile  in  /mnt/'stephenA)ny  must  be  modified  to  use  tlie  correct  bibbny*.f  file 
in  /mnt/stephenA)ny.  These  files  are  used  to  generate  the  elastic  parameiers  and  densities  in 
the  transition  zone  and  are  usually  written  for  each  model.  The  files  created  by  bbny*.bch 


are: 


ROMV*  1  r.  I 


Irva  r\i 


ttrMif  r^r»'*r%rr\y'»>ccrvT* 


BBNY*.LG2  -  log  output  of  bibbny^.f 

BBNY*.LG4  -  log  output  of  the  finite  difference  calculation  including  run  time 
paranreters 

BBNY*.TST  -  binary  file  of  time  scries  at  receiver  locations 
BBN  Y*.SNS  -  binaiy  files  of  the  snapshot  values  (the  last  four  values  before  the 
p>eri(xl  give  the  timestep  value  divided  by  10) 

All  models  use  a  pressure  source  function  which  is  the  third  derivative  of  a 
Gaussian  pulse  with  a  peak  frequency  of  10  Hz.  A  discussion  of  this  pulse  is  given  in 


Appendix  E  of  Stephen  et  al  (1985).  BBNYl  was  run  for  15,000  time  steps  (15  sec.)  and 
all  other  models  were  run  for  7,500  timesteps  (7.5  sec). 

Tlie  layout  ot  tlie  first  model  is  given  in  Figure  1  and  the  velocity -depth  functions 
and  density  depth  functions  arc  given  in  Figure  2  and  Table  1.  The  pressure  lime  series  for 
a  line  of  receivers  at  4.98  km  depth  in  BBNYl  (Figure  3)  show  suable  results  out  to  15  sec 
for  the  primary  bottom  interaction  and  the  first  water  multiple.  At  large  offsets  (greater  than 
3  km)  there  are  shear  wave  peg  leg  multiples  in  the  sediments.  A  weak  reflection  from  the 
absorbing  boundary  can  be  identified  after  13  seconds  on  the  shon  range  traces. 

Figure  4  shows  just  the  primary  bottom  interaction  for  BBNYl.  Tlie  large  first 
arrival  is  tlic  direct  wave  from  the  source.  Tlie  waveform  vaiies  because  of  the  Lloyd's 
mirror  effect  with  the  sea  surface.  A  weak  first  arrival  is  observ'ed  at  ranges  beyond  7.0 
km  which  is  the  head  wave.  About  1.0  second  behind  the  direct  wave  a*  shon  range  is  the 
seafloor  reflection  and  about  0.5  second  behind  this  is  the  basement  reflection.  Tliis  is  all 
Uiat  is  obseived  at  short  range  except  for  some  weak  intra-bed  multiples.  At  larger  ranges 
there  is  later  energy  corresponding  to  shear  wave  multiples  in  the  sediment.  These  data 
could  be  pioccsseJ  fuUhcr  using  Conventional  analysis  techniques  to  bring  out  more  details 
but  this  is  beyond  the  scope  of  this  study, 

ITic  second  model  (BBNY2)  is  the  same  as  the  first  except  the  seafloor  varies 
sinusoidally  with  an  amplitude  of  100  m  and  a  wavelength  of  2.0  km.  Tne  fust  liiil  is 
directly  below  the  shot.  Tlie  parameters  at  the  seafloor  remain  constant  (as  for  BBNYl) 
and  they  vary  linearly  down  to  their  values  just  above  basement.  The  time  series  in 
Figure  5  only  vary  slightly  from  BBNYl  because  they  are  dominated  by  the  direct  wave. 
Tlie  effects  of  bottom  roughness  can  be  seen  in  the  amplitude  variations  of  the  seafloor 
reflection  and  head  wave.  The  shear  wave  multiples  have  lost  their  coherence. 

Tlie  third  model  (BBNY3)  has  a  flat  seaflooi  but  a  sinusoidally  varying  basement 
with  an  amplitude  of  50  m  and  a  wavelength  of  3.0  km.  The  first  valley  is  directly  undei 
the  shot.  Again  values  on  the  interfaces  arc  constant  with  linear  gradients  mnning  vertically 
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.  The  layout  in  two  dimensions  of  the  finite  difference  grid  used  in  this  study. 


Figure  2.  Summary  of  the  velocit  /-depth  and  density-depth  functions  used  for  models  BBNYl 
and  BBNY6. 


Figure  3.  Time  series  resuJts  of  the  finite  difference  code  for  a  line  of  pressure  receivers  at  4.98 
km  depth  in  model  BPNYl.  The  time  series  extend  to  15  seconds  and  show  the  primary  and  first 
water  multiple.  Figure  3a)  has  a  gain  of  1.0  x  10'^  and  shows  the  direct  w.ave  undipped. 

Figure  3b)  has  a  gain  of  1.25  x  10'^  and  shows  more  of  the  low  amplitude  features  such  as  the 
compressional  head  wave  and  the  wide  angle  peg-leg  multiples  of  shear  waves  in  the  sediments. 


model  BBNYl 
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between  the  interfaces.  Small  changes  in  the  amplitude  and  arrival  time  of  the  basement 
reflection  occur  (Figure  6). 

Tlie  fourth  model  (BBNY4)  combines  the  sinusoidal  seafloor  and  basement  of 
BBNY2  and  BBNY3.  The  seafloor  reflection  is  tlie  same  as  for  BBNY2  but  the  basement 
reflection  varies  from  both  BBNY3  and  BBNY4  as  expected  (Figure  7).  'fhe  snapshots  for 
this  relatively  complex  model  should  be  particularly  exciting. 

'rhe  fiftli  model  (BBNY5)  is  a  flat  model,  as  BBNYl,  but  it  has  a  high  velocity 
stringer  Vp  =  5,0  k/s,  Vs  =  2.88  kjs,  and  rho  =  2.0  gm/cc)  between  5.74  and  5.78  km 


depth  and  out  to  a  range  of  0.5  km  (see  Fig.  2).  When  tlris  was  originally  run  with  a  sharp 
contrast  between  the  sediments  and  the  stringer,  the  calculations  were  unstable.  1  do  not 
understand  Uiis  but  it  was  fixed  by  averaging  the  parameters  on  the  boundai-y  values  around 
the  stringer.  The  time  series  on  Figure  8  show  a  large  reflection  from  this  stringer  masking 
basement  at  short  range.  At  a  range  of  0.8  km  the  stringer  is  not  influencing  the  trace.  No 
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unaffected  by  the  stringer. 

ITic  sixth  model  (BBNY6)  dcnxinstrates  the  effect  of  varying  just  the  shear  profile 
in  a  tlat  model  such  as  BBNYl.  Compressional  velocity  and  density  are  unchanged  hut 
shear  velocity  is  increased  making  both  the  sediments  and  basement  more  rigid  (Figure  2). 
In  BBNYl  Poisson's  ratio  in  sediment  and  basement  was  0.46  and  about  0.25, 
respectively.  In  BBNY6  the  coiTcsponding  Poisson's  ratios  are  0.30  and  about  0.05 
(Table  1).  At  short  ranges  the  changes  are  insignificant  but  the  coda  have  larger  amplitude 
at  large  offsets  due  to  wide  angle  shear  wave  reflections  and  multiples  (Irigure  9).  A  nice 
shear  hend  wave  is  also  evident  just  before  the  direct  wave  at  ranges  greater  than  8.4  km. 


Time  series  results  for  a  line  of  pressure  receivers  at  4.98  km  depth  in  model  BBNY3 
seafloor  but  undulating  basement. 


•c  7.  Time  series  n 
undulating  surface 
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■e  8.  Time  series  results  for  a  line  of  pressure  receivers  at  4.98  km  depth  in  model  BBNY5. 
is  a  flat  model  with  a  high  velocity  stringer  just  above  basement.  Reflections  from  the  stringer 
le  seen  at  about  4.5  secs  on  the  0.0  and  0.4  km  traces. 


iglire  9.  Time  series  results  for  a  line  of  pressure  receivers  at  4.98  km  depth  in  mode!  BBNY6 
;ee  Figure  2).  This  mtxlel  is  more  rigid  than  BBNYl  and  shows  significantly  different  results  at 


CONCLUSIONS 


The  results  of  the  test  models  show  that  the  finite  difference  code  at  NORDA  is 
stable  and  producing  reasonable  results  for  large  scale  models  including  the  whole  water 
column.  Further  work  should  involve  applications  of  the  code  to  speeific  problems. 
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